Preliminaries

Loading exec control data

exec_lhq <- as.data.table(fread("subjectdata/axcpt_lhq_data_forpub.csv"))
axcpt <- as.data.table(fread("rawdata/axcpt_trial_data.csv", stringsAsFactors = F, 
    sep = ","))

Create folds

k = 10
k_loo = length(unique(exec_lhq$participant))

exec_lhq <- fold(exec_lhq, k = k, id_col = "participant") %>% ungroup()
exec_lhq <- fold(exec_lhq, k = k_loo, id_col = "participant") %>% ungroup()
## Warning in fold(exec_lhq, k = k_loo, id_col = "participant"): Found 1 existing
## fold column. It will NOT be replaced. Change 'handle_existing_fold_cols' to
## 'remove' if you want to replace it.
exec_lhq <- exec_lhq %>% rename(.folds.10 = .folds_1, .folds.loo = .folds_2)

Entropy PCA and scaling of variables

exec_lhq_ppd <- preprocess_steps(exec_lhq)

Join LHQ and AXCPT

# Merge behavioral and subject data (restricting to only those in subject data
# frame)
axcpt_ppd <- left_join(axcpt, exec_lhq_ppd, by = "participant", suffix = c(".x", 
    ""))
axcpt_ppd <- axcpt_preprocess(axcpt_ppd)

sd.rt <- axcpt_ppd %>% group_by(participant) %>% summarise(meanRT = mean(RT_correct, 
    na.rm = T)) %>% ungroup()

sd.acc <- axcpt_ppd %>% group_by(participant) %>% summarise(meanAcc = mean(accuracy_target))

sd.rt <- sd.rt %>% mutate(gmean = mean(meanRT), gsd = sd(meanRT), lower = gmean - 
    2.5 * gsd, upper = gmean + 2.5 * gsd, outlier = ifelse(meanRT >= lower & meanRT <= 
    upper, "fine", "outlier"))

exclude <- sd.rt$participant[sd.rt$outlier == "outlier"]

axcpt_ppd <- axcpt_ppd[!(axcpt_ppd$participant %in% exclude), ]

exec_lhq <- exec_lhq[!(exec_lhq$participant %in% exclude), ]
exec_lhq_ppd <- exec_lhq_ppd[!(exec_lhq_ppd$participant %in% exclude), ]

exec_lhq_ppd <- fold(exec_lhq_ppd, k = k_loo - length(exclude), id_col = "participant", 
    handle_existing_fold_cols = "remove") %>% ungroup()

k_loo = k_loo - length(exclude)

exec_lhq_ppd <- exec_lhq_ppd %>% rename(.folds.loo = .folds)

newfolds <- exec_lhq_ppd %>% select(participant, .folds.loo)

axcpt_ppd <- axcpt_ppd %>% select(-.folds.loo)
axcpt_ppd <- axcpt_ppd %>% left_join(newfolds)
## Joining, by = "participant"

Jobs data

axcpt_ppd$occupation <- factor(axcpt_ppd$occupation)
contrasts(axcpt_ppd$occupation) <- contr.sum(4)

Experiment version

axcpt_ppd$ExperimentName <- factor(axcpt_ppd$ExperimentName)
contrasts(axcpt_ppd$ExperimentName) <- contr.sum(3)

Data plots

Summary stats table

data.agg <- exec_lhq %>% select(participant, age, aoa, yearsBilingual, L1_current_exposure, 
    L2_current_exposure, L3_current_exposure, ReadEntropy, SpeakEntropy, HomeEntropy, 
    WorkEntropy, SocialEntropy, meanAccent, meanL2ability)

stargazer(data.agg, type = "html", out = "tables/summary_stats.html")
## 
## <table style="text-align:center"><tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Statistic</td><td>N</td><td>Mean</td><td>St. Dev.</td><td>Min</td><td>Pctl(25)</td><td>Pctl(75)</td><td>Max</td></tr>
## <tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">participant</td><td>455</td><td>758.965</td><td>323.781</td><td>212</td><td>453.5</td><td>997</td><td>1,264</td></tr>
## <tr><td style="text-align:left">age</td><td>455</td><td>22.754</td><td>3.633</td><td>18</td><td>20</td><td>24</td><td>35</td></tr>
## <tr><td style="text-align:left">aoa</td><td>455</td><td>6.752</td><td>4.334</td><td>0</td><td>4</td><td>10</td><td>28</td></tr>
## <tr><td style="text-align:left">yearsBilingual</td><td>455</td><td>16.002</td><td>5.340</td><td>1</td><td>13</td><td>19.6</td><td>34</td></tr>
## <tr><td style="text-align:left">L1_current_exposure</td><td>455</td><td>63.690</td><td>20.576</td><td>5</td><td>50</td><td>80</td><td>98</td></tr>
## <tr><td style="text-align:left">L2_current_exposure</td><td>455</td><td>33.890</td><td>20.070</td><td>2</td><td>20</td><td>50</td><td>95</td></tr>
## <tr><td style="text-align:left">L3_current_exposure</td><td>455</td><td>2.420</td><td>5.648</td><td>0</td><td>0</td><td>1</td><td>40</td></tr>
## <tr><td style="text-align:left">ReadEntropy</td><td>455</td><td>0.595</td><td>0.410</td><td>0.000</td><td>0.286</td><td>0.971</td><td>1.571</td></tr>
## <tr><td style="text-align:left">SpeakEntropy</td><td>455</td><td>0.695</td><td>0.417</td><td>0.000</td><td>0.469</td><td>1.000</td><td>1.585</td></tr>
## <tr><td style="text-align:left">HomeEntropy</td><td>455</td><td>0.602</td><td>0.456</td><td>0.000</td><td>0.000</td><td>0.963</td><td>1.585</td></tr>
## <tr><td style="text-align:left">WorkEntropy</td><td>455</td><td>0.756</td><td>0.367</td><td>0</td><td>0.6</td><td>1.0</td><td>2</td></tr>
## <tr><td style="text-align:left">SocialEntropy</td><td>455</td><td>0.944</td><td>0.276</td><td>0.000</td><td>0.811</td><td>1.000</td><td>1.585</td></tr>
## <tr><td style="text-align:left">meanAccent</td><td>455</td><td>3.625</td><td>1.787</td><td>1.000</td><td>2.000</td><td>5.000</td><td>7.000</td></tr>
## <tr><td style="text-align:left">meanL2ability</td><td>455</td><td>7.319</td><td>1.847</td><td>1.000</td><td>6.300</td><td>8.800</td><td>10.000</td></tr>
## <tr><td colspan="8" style="border-bottom: 1px solid black"></td></tr></table>

Group level

sd.rt <- sem(axcpt_ppd, RT_correct, participant, Condition)
sd.acc <- sem(axcpt_ppd, accuracy_target, participant, Condition)

sd.rt$Condition <- factor(sd.rt$Condition, levels = c("AY", "BX", "BY", "AX"))
sd.acc$Condition <- factor(sd.acc$Condition, levels = c("AY", "BX", "BY", "AX"))

p1 <- ggplot(sd.rt, aes(x = Condition, y = mean_RT_correct, ymin = lower, ymax = upper, 
    fill = Condition)) + geom_col() + geom_errorbar(width = 0.2) + ylab("Mean RT (in ms)") + 
    xlab("Trial type") + labs(fill = "Trial type") + theme_minimal(base_size = 12) + 
    scale_fill_brewer(palette = "Set1")


p2 <- ggplot(sd.acc, aes(x = Condition, y = mean_accuracy_target, ymin = lower, ymax = upper, 
    fill = Condition)) + geom_col() + geom_errorbar(width = 0.2) + ylab("Mean accuracy") + 
    xlab("Trial type") + labs(fill = "Trial type") + theme_minimal(base_size = 12) + 
    scale_fill_brewer(palette = "Set1")

(plot <- ggarrange(p1, p2, ncol = 2, labels = c("A", "B"), common.legend = T, legend = "right", 
    align = "v"))
ggsave(plot = plot, file = "figures/mean_acc_rt.png", width = 12, height = 6, units = "cm")

RT: individual differences

rt.sd <- axcpt_ppd %>% group_by(participant, L2_current_exposure, aoa, PCA_General, 
    PCA_Work, Condition) %>% summarise(meanRT = mean(RT_correct, na.rm = T), sdRT = sd(RT_correct, 
    na.rm = T), N = n(), semRT = sdRT/sqrt(N))


rt.sd$Condition <- factor(rt.sd$Condition, levels = c("AY", "BX", "BY", "AX"))

p1 <- ggplot(rt.sd, aes(x = aoa, y = meanRT, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Reaction time (in ms)") + xlab("L2 AoA") + theme_minimal(base_size = 12) + 
    scale_colour_brewer(palette = "Set1")

p2 <- ggplot(rt.sd, aes(x = L2_current_exposure, y = meanRT, colour = Condition)) + 
    geom_point(alpha = 0.2) + geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + 
    xlab("L2 current exposure") + theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p3 <- ggplot(rt.sd, aes(x = PCA_General, y = meanRT, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Reaction time (in ms)") + xlab("PC: General entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p4 <- ggplot(rt.sd, aes(x = PCA_Work, y = meanRT, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Reaction time (in ms)") + xlab("PC: Work entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")


(plot <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, labels = c("A", "B", "C", 
    "D"), common.legend = T, legend = "right", align = "hv"))
ggsave(plot = plot, file = "figures/rt_by_id.png", width = 17.8, height = 17.5, units = "cm")

RT: individual differences: raw

raw <- axcpt_ppd[]
raw$Condition <- factor(raw$Condition, levels = c("AY", "BX", "BY", "AX"))



p1 <- ggplot(raw, aes(x = aoa, y = RT_correct, colour = Condition)) + geom_point(alpha = 0.05) + 
    geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + xlab("L2 AoA") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p2 <- ggplot(raw, aes(x = L2_current_exposure, y = RT_correct, colour = Condition)) + 
    geom_point(alpha = 0.05) + geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + 
    xlab("L2 current exposure") + theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p3 <- ggplot(raw, aes(x = PCA_General, y = RT_correct, colour = Condition)) + geom_point(alpha = 0.05) + 
    geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + xlab("PC: General entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p4 <- ggplot(raw, aes(x = PCA_Work, y = RT_correct, colour = Condition)) + geom_point(alpha = 0.05) + 
    geom_smooth(method = "lm") + ylab("Reaction time (in ms)") + xlab("PC: Work entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")


(plot <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, labels = c("A", "B", "C", 
    "D"), common.legend = T, legend = "right", align = "hv"))
## Warning: Removed 5051 rows containing non-finite values (stat_smooth).
## Warning: Removed 5051 rows containing missing values (geom_point).
## Warning: Removed 5051 rows containing non-finite values (stat_smooth).
## Warning: Removed 5051 rows containing missing values (geom_point).
## Warning: Removed 5051 rows containing non-finite values (stat_smooth).
## Warning: Removed 5051 rows containing missing values (geom_point).
## Warning: Removed 5051 rows containing non-finite values (stat_smooth).
## Warning: Removed 5051 rows containing missing values (geom_point).
## Warning: Removed 5051 rows containing non-finite values (stat_smooth).
## Warning: Removed 5051 rows containing missing values (geom_point).
ggsave(plot = plot, file = "figures/rt_by_id_raw.png", width = 17.8, height = 17.5, 
    units = "cm")

ACC: individual differences

acc.sd <- axcpt_ppd %>% group_by(participant, L2_current_exposure, aoa, PCA_General, 
    PCA_Work, Condition) %>% summarise(meanACC = mean(accuracy_target, na.rm = T), 
    sdACC = sd(accuracy_target, na.rm = T), N = n(), semACC = sdACC/sqrt(N))

p1 <- ggplot(acc.sd, aes(x = aoa, y = meanACC, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Mean accuracy") + xlab("L2 AoA") + theme_minimal(base_size = 12) + 
    scale_colour_brewer(palette = "Set1")

p2 <- ggplot(acc.sd, aes(x = L2_current_exposure, y = meanACC, colour = Condition)) + 
    geom_point(alpha = 0.2) + geom_smooth(method = "lm") + ylab("Mean accuracy") + 
    xlab("L2 current exposure") + theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p3 <- ggplot(acc.sd, aes(x = PCA_General, y = meanACC, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Mean accuracy") + xlab("PC: General entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")

p4 <- ggplot(acc.sd, aes(x = PCA_Work, y = meanACC, colour = Condition)) + geom_smooth(method = "lm") + 
    geom_point(alpha = 0.2) + ylab("Mean accuracy") + xlab("PC: Work entropy") + 
    theme_minimal(base_size = 12) + scale_colour_brewer(palette = "Set1")


(plot <- ggarrange(p1, p2, p3, p4, nrow = 2, ncol = 2, labels = c("A", "B", "C", 
    "D"), common.legend = T, legend = "right", align = "hv"))
ggsave(plot = plot, file = "figures/acc_by_id.png", width = 17.8, height = 17.5, 
    units = "cm")

Set loo fold

axcpt_ppd$.folds <- axcpt_ppd$.folds.loo
axcpt_ppd <- axcpt_ppd[axcpt_ppd$Condition != "AX", ]

Traditional AXCPT analysis

mod.rt.0 <- "RT_correct  ~  Condition + (1|participant)"

mod.rt.1 <- "RT_correct  ~  Condition*aoa + Condition*L2_exposure +  (1|participant)"

mod.rt.1.1 <- "RT_correct  ~  Condition*PCA_Work + Condition*PCA_General +  (1|participant)"

mod.rt.2 <- "RT_correct  ~  Condition*aoa + Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General +  (1|participant)"

mod.rt.3 <- "RT_correct  ~  aoa * (Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General) +  (1|participant)"



mod.acc.0 <- "acc        ~   Condition + (1|participant)"

mod.acc.0.1 <- "acc        ~   Condition + L2_exposure*aoa + (1|participant)"

mod.acc.1 <- "acc        ~   Condition*aoa + Condition*L2_exposure + (1|participant)"

mod.acc.1.1 <- "acc        ~   Condition*PCA_Work + Condition*PCA_General + (1|participant)"

mod.acc.2 <- "acc        ~   Condition*aoa + Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General + (1|participant)"

mod.acc.3 <- "acc        ~   aoa * (Condition*L2_exposure + Condition*PCA_Work + Condition*PCA_General) + (1|participant)"

RT models

# axcpt_ppd <- axcpt_ppd[!(axcpt_ppd$participant %in% exclude),]
axcpt_ppd$Condition <- relevel(axcpt_ppd$Condition, ref = "BX")
all.rt.0 <- lmer(mod.rt.0, axcpt_ppd, REML = F)
all.rt.1 <- lmer(mod.rt.1, axcpt_ppd, REML = F)
all.rt.1.1 <- lmer(mod.rt.1.1, axcpt_ppd, REML = F)
all.rt.2 <- lmer(mod.rt.2, axcpt_ppd, REML = F)
all.rt.3 <- lmer(mod.rt.3, axcpt_ppd, REML = F)
anova(all.rt.0, all.rt.1, all.rt.2, all.rt.3)
## Data: axcpt_ppd
## Models:
## all.rt.0: mod.rt.0
## all.rt.1: mod.rt.1
## all.rt.2: mod.rt.2
## all.rt.3: mod.rt.3
##          Df    AIC    BIC  logLik deviance   Chisq Chi Df Pr(>Chisq)  
## all.rt.0  5 254753 254793 -127372   254743                            
## all.rt.1 11 254754 254841 -127366   254732 11.0123      6    0.08800 .
## all.rt.2 17 254752 254886 -127359   254718 13.8777      6    0.03103 *
## all.rt.3 26 254764 254969 -127356   254712  6.0874      9    0.73114  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
kable(anova(all.rt.0, all.rt.1, all.rt.2, all.rt.3))
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
all.rt.0 5 254753.1 254792.6 -127371.6 254743.1 NA NA NA
all.rt.1 11 254754.1 254840.9 -127366.1 254732.1 11.012288 6 0.0879974
all.rt.2 17 254752.2 254886.4 -127359.1 254718.2 13.877692 6 0.0310327
all.rt.3 26 254764.2 254969.4 -127356.1 254712.2 6.087412 9 0.7311388
kable(anova(all.rt.0, all.rt.1.1, all.rt.2, all.rt.3))
Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
all.rt.0 5 254753.1 254792.6 -127371.6 254743.1 NA NA NA
all.rt.1.1 11 254746.0 254832.9 -127362.0 254724.0 19.095428 6 0.0040055
all.rt.2 17 254752.2 254886.4 -127359.1 254718.2 5.794552 6 0.4465938
all.rt.3 26 254764.2 254969.4 -127356.1 254712.2 6.087412 9 0.7311388
AIC(all.rt.0, all.rt.1, all.rt.1.1, all.rt.2, all.rt.3)
##            df      AIC
## all.rt.0    5 254753.1
## all.rt.1   11 254754.1
## all.rt.1.1 11 254746.0
## all.rt.2   17 254752.2
## all.rt.3   26 254764.2
BIC(all.rt.0, all.rt.1, all.rt.1.1, all.rt.2, all.rt.3)
##            df      BIC
## all.rt.0    5 254792.6
## all.rt.1   11 254840.9
## all.rt.1.1 11 254832.9
## all.rt.2   17 254886.4
## all.rt.3   26 254969.4

So entropy explains more then L2 exp and aoa but not vice versa

summ(all.rt.1.1, confint = T, digits = 3)
Observations 19796
Dependent variable RT_correct
Type Mixed effects linear regression
AIC 254746.033
BIC 254832.859
Pseudo-R² (fixed effects) 0.136
Pseudo-R² (total) 0.315
Fixed Effects
Est. 2.5% 97.5% t val. d.f. p
(Intercept) 384.985 377.229 392.742 97.283 607.866 0.000
ConditionAY 113.831 108.715 118.947 43.608 19364.317 0.000
ConditionBY -42.140 -47.041 -37.240 -16.854 19351.923 0.000
PCA_Work 9.540 1.337 17.743 2.279 608.076 0.023
PCA_General -4.424 -12.618 3.770 -1.058 605.252 0.290
ConditionAY:PCA_Work -5.449 -10.860 -0.037 -1.974 19361.250 0.048
ConditionBY:PCA_Work -1.739 -6.924 3.446 -0.657 19351.850 0.511
ConditionAY:PCA_General 9.271 3.889 14.654 3.376 19363.649 0.001
ConditionBY:PCA_General 7.448 2.288 12.608 2.829 19352.413 0.005
p values calculated using Satterthwaite d.f.
Random Effects
Group Parameter Std. Dev.
participant (Intercept) 74.973
Residual 146.380
Grouping Variables
Group # groups ICC
participant 455 0.208
summary(all.rt.1.1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
##   method [lmerModLmerTest]
## Formula: mod.rt.1.1
##    Data: axcpt_ppd
## 
##       AIC       BIC    logLik  deviance  df.resid 
##  254746.0  254832.9 -127362.0  254724.0     19785 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1265 -0.6608 -0.1892  0.4511  4.8437 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev.
##  participant (Intercept)  5621     74.97  
##  Residual                21427    146.38  
## Number of obs: 19796, groups:  participant, 455
## 
## Fixed effects:
##                          Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)               384.985      3.957   607.866  97.283  < 2e-16 ***
## ConditionAY               113.831      2.610 19364.317  43.608  < 2e-16 ***
## ConditionBY               -42.140      2.500 19351.923 -16.854  < 2e-16 ***
## PCA_Work                    9.540      4.185   608.076   2.279 0.022984 *  
## PCA_General                -4.424      4.181   605.252  -1.058 0.290400    
## ConditionAY:PCA_Work       -5.449      2.761 19361.250  -1.974 0.048452 *  
## ConditionBY:PCA_Work       -1.739      2.646 19351.850  -0.657 0.511042    
## ConditionAY:PCA_General     9.271      2.746 19363.649   3.376 0.000737 ***
## ConditionBY:PCA_General     7.448      2.633 19352.413   2.829 0.004675 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) CndtAY CndtBY PCA_Wr PCA_Gn CAY:PCA_W CBY:PCA_W CAY:PCA_G
## ConditionAY -0.313                                                          
## ConditionBY -0.326  0.493                                                   
## PCA_Work    -0.010  0.013  0.014                                            
## PCA_General  0.005 -0.007 -0.008 -0.331                                     
## CndAY:PCA_W  0.013 -0.041 -0.021 -0.313  0.103                              
## CndBY:PCA_W  0.014 -0.021 -0.042 -0.326  0.107  0.493                       
## CndAY:PCA_G -0.007  0.018  0.011  0.103 -0.311 -0.333    -0.162             
## CndBY:PCA_G -0.008  0.011  0.025  0.107 -0.324 -0.162    -0.331     0.493
anova(all.rt.1.1)
## Type III Analysis of Variance Table with Satterthwaite's method
##                         Sum Sq  Mean Sq NumDF   DenDF   F value    Pr(>F)    
## Condition             82348357 41174178     2 19356.6 1921.5931 < 2.2e-16 ***
## PCA_Work                 72307    72307     1   453.9    3.3746  0.066864 .  
## PCA_General               1873     1873     1   453.0    0.0874  0.767652    
## Condition:PCA_Work       86262    43131     2 19355.2    2.0129  0.133627    
## Condition:PCA_General   282581   141291     2 19356.4    6.5940  0.001372 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ef <- as.data.frame(Effect(c("PCA_Work", "Condition"), all.rt.1.1, xlevels = list(PCA_Work = c(-2, 
    -1, 0, 1, 2)), confidence.level = 0.686))

f1a <- ggplot(ef, aes(x = PCA_Work, y = fit, ymin = lower, ymax = upper, group = Condition)) + 
    geom_line(aes(colour = Condition), size = 1) + geom_ribbon(fill = "grey", alpha = 0.3) + 
    coord_cartesian(ylim = c(300, 550)) + ylab("Model-estimated\nreaction time (in ms)") + 
    xlab("PC: Work Entropy") + scale_colour_brewer(palette = "Set1") + theme_bw(base_size = 12) + 
    labs(color = "Trial type") + theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

# ggsave('figures/axcpt_rt_work_condition.png', width=8.9, height=6, units =
# 'cm')

ef <- as.data.frame(Effect(c("PCA_General", "Condition"), all.rt.1.1, xlevels = list(PCA_General = c(-2, 
    -1, 0, 1, 2)), confidence.level = 0.686))

f1b <- ggplot(ef, aes(x = PCA_General, y = fit, ymin = lower, ymax = upper, group = Condition)) + 
    geom_line(aes(colour = Condition), size = 1) + geom_ribbon(fill = "grey", alpha = 0.3) + 
    coord_cartesian(ylim = c(300, 550)) + ylab("Model-estimates\nreaction time (in ms)") + 
    xlab("PC: General Entropy") + labs(color = "Trial type") + scale_colour_brewer(palette = "Set1") + 
    theme_bw(base_size = 12) + theme(panel.border = element_blank(), panel.grid.major = element_blank(), 
    panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

# ggsave('figures/axcpt_rt_general_condition.png', width=8.9, height=6, units =
# 'cm')

(f1 <- ggarrange(f1b, f1a, labels = c("A", "B"), legend = "bottom", common.legend = T))

ggsave(plot = f1, "figures/traditional_rt_lme4_effects.png", width = 17.8, height = 8, 
    units = "cm")

RT models cross validated

cv.all.rt.0 <- crossvalidate(mod = mod.rt.0, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)
## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"
cv.all.rt.1 <- crossvalidate(mod = mod.rt.1, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)
## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"
cv.all.rt.1.1 <- crossvalidate(mod = mod.rt.1.1, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)
## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"
cv.all.rt.2 <- crossvalidate(mod = mod.rt.2, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)
## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"
cv.all.rt.3 <- crossvalidate(mod = mod.rt.3, k = k_loo, data = axcpt_ppd, dependent = "RT_correct", 
    dv_continuous = T, random = TRUE, returnRuns = T)
## [1] "current iteration: 10 / 455"
## [1] "current iteration: 20 / 455"
## [1] "current iteration: 30 / 455"
## [1] "current iteration: 40 / 455"
## [1] "current iteration: 50 / 455"
## [1] "current iteration: 60 / 455"
## [1] "current iteration: 70 / 455"
## [1] "current iteration: 80 / 455"
## [1] "current iteration: 90 / 455"
## [1] "current iteration: 100 / 455"
## [1] "current iteration: 110 / 455"
## [1] "current iteration: 120 / 455"
## [1] "current iteration: 130 / 455"
## [1] "current iteration: 140 / 455"
## [1] "current iteration: 150 / 455"
## [1] "current iteration: 160 / 455"
## [1] "current iteration: 170 / 455"
## [1] "current iteration: 180 / 455"
## [1] "current iteration: 190 / 455"
## [1] "current iteration: 200 / 455"
## [1] "current iteration: 210 / 455"
## [1] "current iteration: 220 / 455"
## [1] "current iteration: 230 / 455"
## [1] "current iteration: 240 / 455"
## [1] "current iteration: 250 / 455"
## [1] "current iteration: 260 / 455"
## [1] "current iteration: 270 / 455"
## [1] "current iteration: 280 / 455"
## [1] "current iteration: 290 / 455"
## [1] "current iteration: 300 / 455"
## [1] "current iteration: 310 / 455"
## [1] "current iteration: 320 / 455"
## [1] "current iteration: 330 / 455"
## [1] "current iteration: 340 / 455"
## [1] "current iteration: 350 / 455"
## [1] "current iteration: 360 / 455"
## [1] "current iteration: 370 / 455"
## [1] "current iteration: 380 / 455"
## [1] "current iteration: 390 / 455"
## [1] "current iteration: 400 / 455"
## [1] "current iteration: 410 / 455"
## [1] "current iteration: 420 / 455"
## [1] "current iteration: 430 / 455"
## [1] "current iteration: 440 / 455"
## [1] "current iteration: 450 / 455"
cv.rt.0.iter = cv.all.rt.0$performances
cv.rt.0 = cv.all.rt.0$performances_mean

cv.rt.1.iter = cv.all.rt.1$performances
cv.rt.1 = cv.all.rt.1$performances_mean

cv.rt.1.1.iter = cv.all.rt.1.1$performances
cv.rt.1.1 = cv.all.rt.1.1$performances_mean

cv.rt.2.iter = cv.all.rt.2$performances
cv.rt.2 = cv.all.rt.2$performances_mean

cv.rt.3.iter = cv.all.rt.3$performances
cv.rt.3 = cv.all.rt.3$performances_mean
xval_axcpt_rt_k_loo.iter <- rbind.data.frame(cv.rt.0.iter, cv.rt.1.iter, cv.rt.1.1.iter, 
    cv.rt.2.iter, cv.rt.3.iter)
write.csv(xval_axcpt_rt_k_loo.iter, "tables/xval_axcpt_rt_k_loo.iter.csv", row.names = F)

xval_axcpt_rt_k_loo.iter <- xval_axcpt_rt_k_loo.iter %>% group_by(mod) %>% mutate(meanRMSE_test = mean(RMSE_test)) %>% 
    ungroup() %>% mutate(mod = fct_reorder(mod, desc(meanRMSE_test)))

mean.xval_axcpt_rt_k_loo.iter <- xval_axcpt_rt_k_loo.iter %>% group_by(mod) %>% summarise(meanRMSE_train = mean(RMSE_train), 
    meanRMSE_test = mean(RMSE_test), meanr2_train = mean(r2_train), meanr2_test = mean(r2_test), 
    sd = sd(RMSE_test), N = n(), sem = sd/sqrt(N)) %>% ungroup() %>% mutate(mod = fct_reorder(mod, 
    desc(meanRMSE_test)))

write.csv(mean.xval_axcpt_rt_k_loo.iter, "tables/xval_axcpt_rt_k_loo.mean.csv", row.names = F)

ACC models

all.acc.0 <- glmer(mod.acc.0, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.0.1 <- glmer(mod.acc.0.1, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.1 <- glmer(mod.acc.1, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.1.1 <- glmer(mod.acc.1.1, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.2 <- glmer(mod.acc.2, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))

all.acc.3 <- glmer(mod.acc.3, axcpt_ppd, family = "binomial", control = glmerControl(optimizer = "bobyqa", 
    optCtrl = list(maxfun = 3e+05)))
kable(anova(all.acc.0, all.acc.1, all.acc.2, all.acc.3))
Df AIC BIC logLik de viance Chisq Ch i Df Pr (>Chisq)
all.acc.0 4 13487.86 13519.91 -6739.931 13479.86 NA NA NA
all.acc.1 10 13490.65 13570.76 -6735.323 13470.65 9.215386 6 0.1618223
all.acc.2 16 13496.08 13624.26 -6732.038 13464.08 6.570276 6 0.3624199
all.acc.3 25 13495.63 13695.92 -6722.816 13445.63 18.444541 9 0.0303527
kable(anova(all.acc.0, all.acc.1.1, all.acc.2, all.acc.3))
Df AIC BIC logLik de viance Chisq Ch i Df Pr (>Chisq)
all.acc.0 4 13487.86 13519.91 -6739.931 13479.86 NA NA NA
all.acc.1.1 10 13494.19 13574.31 -6737.094 13474.19 5.673929 6 0.4606892
all.acc.2 16 13496.08 13624.26 -6732.038 13464.08 10.111734 6 0.1200248
all.acc.3 25 13495.63 13695.92 -6722.816 13445.63 18.444541 9 0.0303527
kable(anova(all.acc.0, all.acc.0.1))
Df AIC BIC logLik de viance Chisq Ch i Df Pr (>Chisq)
all.acc.0 4 13487.86 13519.91 -6739.931 13479.86 NA NA NA
all.acc.0.1 7 13483.66 13539.74 -6734.830 13469.66 10.20268 3 0.0169196
summ(all.acc.0.1, confint = T, digits = 3)
Observations 22286
Dependent variable acc
Type Mixed effects generalized linear model
Family binomial
Link logit
AIC 13483.659
BIC 13539.741
Pseudo-R² (fixed effects) 0.299
Pseudo-R² (total) 0.345
Fixed Effects
Est. 2.5% 97.5% z val. p
(Intercept) 1.994 1.911 2.077 46.969 0.000
ConditionAY -0.332 -0.423 -0.242 -7.162 0.000
ConditionBY 2.430 2.219 2.641 22.597 0.000
L2_exposure -0.045 -0.109 0.020 -1.353 0.176
aoa 0.015 -0.051 0.081 0.442 0.659
L2_exposure:aoa 0.089 0.027 0.150 2.816 0.005
Random Effects
Group Parameter Std. Dev.
participant (Intercept) 0.484
Grouping Variables
Group # groups ICC
participant 455 0.067
summ(all.acc.0, confint = T, digits = 3)
Observations 22286
Dependent variable acc
Type Mixed effects generalized linear model
Family binomial
Link logit
AIC 13487.862
BIC 13519.909
Pseudo-R² (fixed effects) 0.296
Pseudo-R² (total) 0.345
Fixed Effects
Est. 2.5% 97.5% z val. p
(Intercept) 2.002 1.918 2.085 46.884 0.000
ConditionAY -0.333 -0.424 -0.242 -7.171 0.000
ConditionBY 2.430 2.219 2.640 22.616 0.000
Random Effects
Group Parameter Std. Dev.
participant (Intercept) 0.495
Grouping Variables
Group # groups ICC
participant 455 0.069
anova(all.acc.0)
## Analysis of Variance Table
##           Df Sum Sq Mean Sq F value
## Condition  2 643.32  321.66  321.66
ef <- as.data.frame(Effect(c("L2_exposure", "aoa"), all.acc.0.1, xlevels = list(aoa = c(-2, 
    -1, 0, 1, 2)), confidence.level = 0.686))
ef$aoa <- factor(ef$aoa)
f_acc <- ggplot(ef, aes(x = L2_exposure, y = fit, ymin = lower, ymax = upper, group = aoa)) + 
    geom_line(aes(colour = aoa), size = 1) + geom_ribbon(fill = "grey", alpha = 0.3) + 
    coord_cartesian(ylim = c(0.5, 1)) + ylab("Model-estimated\naccuracy") + xlab("Scaled L2 exposure") + 
    scale_colour_brewer(palette = "Dark2") + theme_bw(base_size = 12) + theme(panel.border = element_blank(), 
    panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"))

ggsave("figures/traditional_acc_plot.png", plot = f_acc, width = 12, height = 8, 
    units = "cm")